Analysis of a spatial Lotka-Volterra model with a finite range predator-prey 

interaction 
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We perform an analysis of a recent spatial version of the classical Lotka-Volterra model, where a 
finite scale controls individuals' interaction. We study the behavior of the predator-prey dynamics 
in physical spaces higher than one, showing how spatial patterns can emerge for some values of the 
interaction range and of the diffusion parameter. 
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I. INTRODUCTION 



Reaction-diffusion theory plays a very important role 
in the study of pattern formation in biology, an aspect of 
populations dynamics which has gained importance since 
some Moran's seminal works [1,2]. More recently, theo- 
retical ecology studies focused on the dynamical mech- 
anisms which can be considered responsible for produc- 
ing such structures. Their principal aim is to determine 
the local rules which cause the emergence of macroscopic 
patterns [7]. At present, a central issue is to determine 
whether the predator-prey interactions can be considered 
among these mechanisms. Actually, spatial correlations 
between predators and preys populations have been ob- 
served in real systems such as predatory beetles and lar- 
val flics as their preys [3] . This research shows that these 
larval flies appear to gain local clustering as beetle abun- 
dance approaches a carrying capacity with the prey pop- 
ulation. 

A traditional approach for describing these systems is 
based on spatial Lotka-Volterra type models. This de- 
scription is able to generate spatial structures of the Tur- 
ing type, where the instabilities are directly driven by 
mechanisms connected with the process of diffusion [10]. 
These works are inspired to the classical ideas introduced 
by Turing, who showed that diffusion can destabilize the 
homogeneous equilibrium solution of a reaction-diffusion 
system [9]. This topic has recently received some atten- 
tion, in particular in relation to the effects of stochas- 
ticity. For example, in [4] it is shown that Turing-like 
patterns can indeed emerge beyond the parameter re- 
gion predicted by the conventional Turing theory. This 
phenomenon is due to the effect of finite size corrections 
on the mean-field idealized dynamics. Also demographic 
noise can induce persistent spatial pattern formation and 
temporal oscillations [5] . Specifically, demographic noise 
greatly enlarges the region of parameters space where 
pattern formation occurs, in relation to the prediction 
of the mean- field theory. Finally, Scott et. al. [6] intro- 
duce an analytic method to describe stochastic spatially 
varying systems governed by reaction-transport master 



equation. 

More recently, a different mechanism capable of gener- 
ating spatial patterns has been explored [13]. It is based 
on the introduction of a probability of interaction (preda- 
tion) which becomes a function of the distance between 
individuals. This description presents two interesting 
points. First, it is not subject to the strong constraints 
which limit the emergence of spatial inhomogeneities in 
the traditional approach [14]. Second, it is based on a 
simple ingredient which has a deep ecological motivation. 
In fact, it has became evident that predation strongly de- 
pend on probability of encounter [8] which must take into 
account predators shifts in response to prey movements. 
Moreover, territoriality is another fundamental aspect in 
the ecology of many predatory animals such as wolves, 
lions, hyenas, African wild dogs and badgers. In this case 
an immediate question arises as to how the predation is 
strongly influenced by a land which is divided up into 
predator territories [10]. These elements result in a spa- 
tial and temporal variation in the predation risk [12]. A 
possible description of one aspect of this spatial depen- 
dence can be obtained considering that the probability 
at which a consumer meets a prey should be dependent 
on the relative distance between them [15]. The intro- 
duction of a spatial scale of interaction has been widely 
applied to model competition in community ecology [15] 
and more recently it was also used in studies of evolu- 
tionary theory [16]. 

In the following, we present a detailed analysis of this 
approach for an implementation of a predator-prey sys- 
tem in a two dimensional physical space, the most rele- 
vant for describing real ecosystems. 



II. THE MEAN-FIELD MODEL 

We start with the mean-field version of the model. 
This description is characterized by a couple of equa- 
tions, one for the preys N(x, t) and other for the preda- 
tors P(x,t). They describe diffusion in real space and 
the strength of the interaction in the nonlinear term is 
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a function of the density of individuals in the proximity 
[17]. The reason for introducing a spatial scale for the 
interaction can be clearer understood considering that, 
based on the reported ecological motivations, we want 
to describe a probability of reproduction/death not con- 
stant but dependent on the number of predators / preys in 
the surrounding. In fact, in a region more populated by 
preys, predators should have a better chance of reproduc- 
tion. Conversely, in a region more populated by preda- 
tors, preys should have a greater chance of dying. This 
nonlocality appears as an integral term which takes into 
account the local density of population. It is important 
to point out that this effect is not related to the individ- 
uals' velocity, controlled by the diffusion coefficients, but 
to the memory which individuals maintain in relation to 
a specific territory. Territoriality is a fundamental aspect 
in the ecology of many predatory animals [10]. For this 
reason, as the range of interaction is not simply linked to 
the relative motion between predators and preys, it does 
not depend on a single interaction scale. 

Considering a two dimensional physical space the equa- 
tions which describe our system are the following: 
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where x = (x, y) is the position of predators or preys. 
Predators consume the preys with an intrinsic rate a and 
reproduce with rate /3, r is the prey's growth rate and 
predators are assumed to spontaneously die with rate m. 
Dn and Dp are the diffusion coefficients of preys and 
predators, respectively. 

A microscopic version of these equations is introduced 
in the next section. There, the integral terms correspond 
to counting the number of preys/predators which are at 
a shorter distance than Ri from the considered preda- 
tor/prey. 

We would like to point out that a rigorous derivation, 
which would show that the continuum field equations in- 
troduced here are totally analogous to the discrete model 
presented in the next section, can be obtained by using 
Fock space techniques after making some usual assump- 
tion and approximation [17,18]. 

To perform this passage between discrete model and 
mean field description the microscopic model should be 
defined on a lattice. The lattice model will be equivalent 
to the off-lattice one in the limit of small lattice size. For 
this reason the off-lattice model present in the next sec- 
tion should be modified to fit into a grid. If we let the grid 
spacing to become negligibly small we recover the contin- 
uum limit. With the aim of preserving the non- locality 



of our macroscopic model, in which the interaction range 
is finite, this operation of limit should be made carefully. 
In fact, only the grid size, but not the interaction lengths 
Ri should vanish. This is obtained fixing Ri to a finite 
value when going to the continuum. 

Analogous models, but not characterized by a finite 
scale of interaction, appeared recently in the literature. 
Spatio-temporal patterns can emerge sustained by the 
stochastic component of the dynamics. In this case, the 
demographic fluctuations, directly reproduced in the in- 
dividual based description level, are the responsible for 
pattern formation [5,19]. 

The one dimensional case was deeply analyzed in the 
reference [13]. In the following we will present some re- 
sults related to the two-dimensional physical space, and, 
for theoretical interest, to n-dimcnsional spaces. 

Also in the two dimensional space, an absorbing phase 
N(x,t) = P(x,t) — and a survival phase N(x,t) = 
n ™ R -i ; P(x,t) ~ nap' 2 ex ist and they correspond to spa- 
tially homogeneous, stationary solutions. 

The analysis for detecting the presence of solutions 
with spatial structure is based on the idea of perturb- 
ing the survival phase stationary solution N(x,t) and 
P{x, t). This is obtained introducing small harmonic per- 
turbations [10,17]: N(x,t) = N + A N exp[Xt + ik ■ x]; 
P(x,t) = P + Apexp[Xt + ik ■ x], where k = (k x ,k y ) 
is a two-dimensional wave vector of modulus \k\ = k. 
The simple harmonic form of these perturbations reduces 
equations 1 to a linear system. The solution of this sys- 
tem generates some constraints on the values which A and 
k can assume. This dispersion relation can be written in 
the following form: 



A(fc) = - — (D N +D P ) 
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where J\ is the first-order Bessel function. 

If, for some k, Re[X(k)] > spatial patterns can 
emerge. It follows that R\ 7^ R2 is a necessary condi- 
tion. Moreover, because for = Dp = the condition 
is always satisfied, it is evident that the instability is 
driven by the range of the interaction and is independent 
of the diffusion process. 

If we consider the simpler situation with Dm = Dp = 
D and i?i = 2R 2 = 2R, and we introduce the rescaled 

variables K = kR and A = A^-, for /j, = y/rmR 2 / D 
equation 2 reduces to: 



\{K) = -K 2 + j- yj -Ji{K)J!(2K). 



(3) 



The edge of spatial patterns emergence corresponds to 
the values of the parameters for which the maximum K m 
of X(K) becomes positive. For the original variables this 
happens if: 



D 



> 22.4228 and 



2.19535 



R 



(4) 



Figure 1 shows X(K) for different ji values. 
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FIG. 1. Dispersion relation X(K) for different values of the 
parameters (p = sJrraB 2 / D). 

In the following we turn our attention to the model 
implemented in a n-dimensional space, with n > 2. In 
this case, the corresponding equations are obviously ob- 
tained generalizing the Laplacian terms and the integra- 
tion domains. Performing a calculation analogous to the 
2-dimcnsional case, we arrive at the following dispersion 
relation: 



X(k) = - — {D N + D P ) 
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where n is the space dimension. Introducing the vari- 
ables K = fci?2, 
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patterns can emerge if the following condition is satisfied: 



p 2 K 4 
p? 



+ 2 n T 1 



r^ 2 J n /2(K)J n/2 (pK) 
p n / 2 K n 



< 0. 



(6) 



As for the two dimensional situation, a critical value p c 
exists for which, if p > p c spatial patterns emerge. In 
Figure 2 we show the p~ x values as a function of the ra- 
tio p = for different dimensions. The value of p\. x 
increases as the parameter p increases. This means that 
a stronger difference between Ri and i?2 leads to easier 
cluster formation. Moreover, it is clear that for p > 1 
clustered distributions can appear also for a physical di- 
mension greater than 2. In higher dimensions cluster 
formation becomes more difficult but it is viable. It 
is important to remember that even the case of a 3- 
dimensional space is not a common ecosystem. Obviously 
it does not exist for terrestrial individuals and even ma- 
rine species tend to remain within a layer of small thick- 
ness compared to their horizontal movements. Anyway, 
our results for dimension higher than 2 have a theoret- 
ical interest. In fact, qualitative differences depending 
on the space dimension are a common feature in systems 
with diffusive phenomena and clusters formation (see for 
example [20]). 
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FIG. 2. p c 1 values in dependence of p for different physical 
dimensions n. 



III. THE MICROSCOPIC DISCRETE MODEL 

In the following we introduce a microscopic discrete 
stochastic formulation of the mean field model. The im- 
portance of comparing the results of the mean-field de- 
scription with an individual based one is due to the pos- 
sibility of introducing demographic fluctuations. These 
fluctuations result in the appearance of an intrinsic 
stochasticity which, in principle, is always present in real 
phenomena. Moreover, this technique can describe the 
occurrence of threshold effects, generated by the discrete 
nature of individuals. These effects are not present in the 
mean field description where every small amount of the 
density of population is acceptable, even if unrealistically 
small [21]. 

We implement our individual based model on the 
square [0, 1] x [0, 1], with periodic boundary conditions. 
Simulations start with an initial population of Pq preda- 
tors and No preys, randomly located. A time step of our 
simulation corresponds to carrying out the following pro- 
cesses: 1) diffusion, where a predator, or a prey, is ran- 
domly selected and moves some distance, in a random 
direction, chosen from a Gaussian distribution of stan- 
dard deviation a. 2) predator reproduction, with rate 
/3-/Vj| , where iVg is the number of preys which are at a 
shorter distance than i?2 from the predator at position 
x. 3) predator death, with probability m. 4) prey re- 
production, with probability r. 5) prey death, with rate 
aP^ , where Pjj is the number of predators which are at 
a shorter distance than R\ from the prey at position y. 
All the newborns maintain the same location as the par- 
ents. We evaluate N^ 2 and Pjj using periodic boundary 
conditions. These processes are executed sequentially by 
the whole population. In fact, for each action, an indi- 
vidual of each population (predator or prey) is randomly 
selected and this operation is repeated for a number of 
times equal to the size of the corresponding population. 

It is important to point out that if we measure time 
in units of the simulation time step, the coefficient D of 
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equation 1 is related to the discrete model through the 
relation D = a 2 /2. Birth and death probabilities are the 
same in the continuous and in the discrete model. 

Analyzing the data generated throughout these simu- 
lations we are able to explore the temporal and spacial 
behavior of our system. The temporal dynamics of our 
simulations is strongly dependent on the form of the spa- 
tial organization of the populations. The temporal evo- 
lution of the simulations characterized by homogeneous 
spatial distributions can be easily understood consider- 
ing a classical Lotka-Volterra model. In fact, if we start 
with an initial population of Pq = predators and 

Ao = ^pfti preys, corresponding to the two stationary 
solutions of the mean field model, inevitable fluctuations 
tend to push the system away from the trivial survival 
phase and induce irregular population oscillations that 
almost resemble the deterministic cycles of the classi- 
cal Lotka-Volterra model. This behaviour has been de- 
scribed also in other analogous stochastic models [22,18]. 
We must remember that an important outcome of the 
introduction of intrinsic stochasticity is the possibility 
of predators' extinction. This fact, combined with the 
possibility of extremely rapid increase in the number of 
preys, force to a careful choice of the initial conditions 
and parameters of our simulational runs. This is nec- 
essary to obtain controlled irregular oscillations, which 
swing in a rather erratic fashion around an average value 
(see Figure 3). Otherwise, the explosion or extinction 
of population do not allow to carry out sufficiently long 
standing simulations. 

The situation changes when the system is evolving in 
the spatial clustered regime. In fact, for some values of 
the parameters, the survival phase of the homogeneous 
distribution is different from the clustered one. For exam- 
ple, starting with Pq = ^^ R -i predators and No = v ™ R s 
preys, after a rapid transient, when spatial structures 
emerge, the oscillations organise around a new survival 
phase. In fact, the difference between P Ri and N R2 be- 
comes irrelevant and the values of a and /3 determine 
the ratio in the new population equilibrium (see Fig- 
ure 4). This behaviour can be easily understood con- 
sidering that the dimension of the bulk of the clusters is 
generally smaller than i?,, as shown at the end of this 
section. We can also remark that the more the spatial 
solution is marked by clustering (lower R and D values), 
the more the irregularities of the temporal oscillations 
increase. It is interesting to point out a recent model 
which shows analogous time oscillations in predators and 
preys number associated with a spatial structure. These 
oscillations are due to stochastic fluctuations about the 
time- independent solutions of the deterministic equations 
[11]. Other models describe similar oscillating behaviour 
[24]. Interestingly, these works identify and analytically 
describe a mechanism of amplification of demographic 
noise which can give rise to coherent oscillations. In par- 
ticular, as in our model, in references [5,19] these oscilla- 
tions are spatiotemporal in nature. 
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FIG. 3. Temporal evolution for populations with homo- 
geneous spatial distributions. Top: Ri = 2Ki = 0.16, 
a = 0.03, r = m = 0.5, a = 6.22 x 10 _4 ,/3 = 1.55 x 10~ 4 , 
P = N = 40000. Bottom: R t = R 2 = 0.08, o = 0.008, 
r = m = 0.5, a = p = 6.22 x 10~ 4 , P = N = 40000. The 
two temporal evolutions are very similar. This analogous be- 
haviour is caused by the fact that no patterns are generated 
in either of the two simulations: in the upper one because of 
the large value of a, in the lower one because of RI = R2. 
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FIG. 4. Top: temporal evolution for populations with 
modulated spatial distributions (Ri = 2R2 = 0.16). 
Other parameters are: a = 0.008, r — m = 0.5, 
a = 6.22 x 10" 4 ,^ = 1.55 x 10" 4 , Po = iVo = 40000. Bot- 
tom: on the left, modulated spatial distribution correspond- 
ing to the simulation shown on the top of the figure. On the 
right, homogeneous distribution generated using the follow- 
ing parameters: Ri — R2 — 0.08, a — 0.008, r = m = 0.5, 
a = /3 = 6.22 x 10~ 4 , P = N = 40000. 

In the following we present some results related to 
the spatial distributions generated by our model. As 
reported in Figure 4 regular patterns can be obtained 
only for simulations where the two ranges of interac- 
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tion, Ri and R2, have a different vaiue. This out- 
come is in agreement with the results obtained by the 
mean field description. The spatial patterns are char- 
acterised by a sequence of isolated colonies, generally 
called of spikes [17,23]. These spikes keep changing along 
the time, according to the oscillation reported in the 
temporal evolution of the total population (top of Fig- 
ure 4). These dynamics result in periodic stationary spa- 
tial structures characterised by fluctuating clusters ar- 
ranged on an hexagonal lattice. We can remark that 
changing the geometry of the domains of integration it is 
possible to generate different clusters arrangements (for 
example, a square domain generates clusters arranged on 
a square lattice). 

In the following, we characterise the transition between 
the homogeneous spatial distribution and the inhomoge- 
neous one (segregation transition). A proper order pa- 
rameter is obtained evaluating a structure factor S(K) 
defined as: 



S(K) 



N{r) 

£ 



1 



■ exp[i<7 



A" 



(7) 



where the sum is carried out over all predators (or preys) 
N(t) at time r, Xj(r) = (xj(r), j/j(t)) is the position 
of the predator j, q = (q x , q y ) is a two-dimensional wave 
vector and the average is a spherical average over all wave 
vectors of modulus |g| = K. The position of the sec- 
ondary local maximum of this function {Km) identifies 
the emergence of periodic structures in the spatial dis- 
tribution. In fact, the transition from a homogeneous to 
an inhomogeneous distribution matches the jump of Km 
to higher values, corresponding to the wave number of 
the periodic clusters present in the space. This is clearly 
shown in the insets of Figure 5. The segregation tran- 
sition is characterised by the passage of Km from small 
values to higher values as soon as a modulation becomes 
dominant [17,23]. 

In Figure 5 we show Km as a function of a. As pre- 
dicted by the mean field approximation, for any value of 
the range of the interaction, a critical value of the diffu- 
sion coefficient exists above which no spatial structures 
emerge. In the figure we show some results obtained fix- 
ing the parameters R\ = 2i?2 = 0.16. For the sake of 
simplicity we explored the case with i?i = 2i?2, but sim- 
ilar behaviour can be obtained for other values of p 7^ 1 . 

A similar analysis can be unfolded fixing a and looking 
for the R dependence. We obtained spatial modulations 
of arbitrary wavelengths just tuning the parameter R2. 
In Figure 6 we can see how the continuous description 
can reproduce quantitatively the period of the patterns 
which emerges from the modulated distributions of the 
Monte Carlo simulations. In fact, the wavenumber of 
these distributions are well approximated by the second 
relation in equation 4, which gives the fastest growing 
mode of the mean-field approximation. For extremely 
short-ranged interaction, a noisy spatially homogeneous 
distribution appears. This means that clusters appear for 



i?2 > R c , where R c is a critical value of the interaction 
length for which the segregation transition takes place. 

It is interesting to remark another phenomenon. Also 
in a two dimension implementation, independently of the 
values of i?, for low D values and low population density 
the demographic fluctuations can introduce a source of 
spatial correlation which can generate disordered clusters 
[13]. 
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FIG. 5. K M as a function of a (R x = 2R 2 = 0.16, 
r = m = 0.5, a = 6.22 x 10 _4 ,/3 = 1.55 x 10" 4 , 
Po = No = 40000). For this parameters the analytical pre- 
diction gives a c = 0.019, in good accordance with the data 
obtained from the discrete model. 
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FIG. 6. Km as a function of R 2 {a = 0.004, 
r = m = 0.5,.Ri = 2R 2 , P = N = 80000, a = [27rP -R§]"\ 
j3 = 0.25a). Simulational data are compared with the ana- 
lytical predictions given by eq. 4, from which we obtained: 
k m ss 2.2/i?2 (continuous line) and R c R3 0.019 (dashed line). 

Finally, we try to evaluate the typical size S of the clus- 
ters which appear in the spiky phase. The typical cluster 
size was calculated averaging the root mean square dis- 
persion for all the clusters present in the system at a given 
time. Data exposed in Figure 7 show a dependence of the 
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cluster size on the diffusion coefficient: S oc er, equivalent 
to S oc \/~D. We can easily interpret this result. Since 
there is no attraction or repulsion between individuals, 
they will experience just a diffusive motion during their 
lifetimes. Assuming that the individuals confined in a 
cluster diffuse a distance proportional to y/D, the de- 
pendence of the clusters size on the diffusion coefficient 
shown in Figure 7 is straightforward. Similar results were 
reported for simpler models in [23,25]. 



| 0.03- 




0025 005 0.0075 0.0 ] 00125 

a 

FIG. 7. Cluster size as a function of a. The dashed line 
has slope 1, which corresponds to a square root dependence 
of the cluster size on the diffusion coefficient. 



IV. CONCLUSIONS 

We have developed an in-deep analysis of a spatial 
Lotka-Volterra model characterised by a finite range 
predator-prey interaction in a two dimensional space. 
This physical space is obviously the most appropriate for 
describing real ecosystems. 

This model can be considered relevant from an eco- 
logical perspective for two reasons. First, because it con- 
tributes to the discussion about the mechanisms that can 
generate spatial correlations in population ecology, show- 
ing how preys and predators interaction can be consid- 
ered among them. This conjecture is of particular im- 
portance as some records of spatial correlations between 
preys and predators were clearly reported [2,3]. Second, 
because the emergence of spatial patterns in an erratic 
oscillatory regime, which can contemplate predators ex- 
tinction, shows realistic elements generally absent from 
conventional descriptions. 

A theoretical interest related to this kind of studies 
exists too. From this perspective, we demonstrate that 
in this model instability is driven by the form of the 
interaction and is independent of the diffusion process. 
Moreover, we carried out our analysis developing a direct 
comparison between an individual based implementation 
and a mean-field description, showing a perfect match 
between them for many quantitative features, a fact that 



is not always achieved [5,11,23,26]. Finally, we show how 
the presence of structures persists independently of the 
space dimension where the model can be implemented. 
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